There are three basic steps for constructing plots with ggplot2:

1. Supply a data set you want to plot to ggplot().
2. Define an aesthetic mapping with aes().
   This describes how variables of your data are mapped to visual properties, e.g. variable "age" is plotted on the x-axis and "income" on the y-axis.
3. Add layers of geoms (geometrical objects) that describe which graphical primitives to use (e.g. points in a scatter plot or bars in a bar plot).

Additionally, you can further change the appearance of your plot by:

• altering the scales (e.g. use a logarithmic scale, modify display of factors, etc.)
• defining facets → create small multiples, each showing a different subset of the data
• changing the overall appearance of the plot by adjusting its theme (e.g. change background color, rotate axis labels, etc.)

You combine all these steps with a +.

We usually work with two-dimensional vector data: Points located within a coordinate reference system (CRS).

• the point -73.94°, 40.6° is represented in the WGS84 CRS: most popular world coord. system (e.g. used by GPS)
• this CRS is ellipsoidal; its components are given in degrees:
  • -73.94° is the longitude: about 74° West from the meridian
  • 40.6° is the latitude: about 41° North from the equator

Let's put this point into context:

• points can be connected to form other geometrical shapes such as lines or polygons
• points and shapes can represent numerous geographic entities: cities, buildings, countries, rivers, lakes, etc.

Most packages for working with geo-data in R rely on the Simple Features standard (→ sf package).

Simple features describe how objects in the real world can be represented in computers in terms of their spatial geometry:

Features have a geometry describing where on Earth the feature is located, and they have attributes, which describe other properties. – sf package vignette

Examples:

• geometry: point at -73.94°, 40.6°
• name: New York City
• population: 8.6 mio.

or

• geometry: polygons consisting of points at …
• name: Italy
• most popular meal: pizza

A spatial dataset is a dataset that contains geo-data (in a geometry column) and additional attributes of that entity like population, poverty rate, etc.

First, we load the packages that we need:

library(ggplot2)
library(maps)
library(sf)

The function map can be used to load the "world" data. We need to convert it to a "simple features" object via st_as_sf():

worldmap_data <- st_as_sf(map('world', plot = FALSE, fill = TRUE))
head(worldmap_data, n = 3)  # to have a look inside

## Simple feature collection with 3 features and 1 field
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -8.68335 ymin: 18.98662 xmax: 74.89131 ymax: 42.64795
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
##                         geometry          ID
## 1 MULTIPOLYGON (((74.89131 37... Afghanistan
## 2 MULTIPOLYGON (((20.06396 42...     Albania
## 3 MULTIPOLYGON (((8.576563 36...     Algeria

• dataset with a "header" giving some general information
• two columns: geometry (spatial data) and ID ("meta data")

What do the numbers in the geometry column mean? What's their CRS?

We have some cities along with their longitude (lng) and latitude (lat):

##       name    lng    lat
## 1   Berlin  13.38  52.52
## 2 New York -73.94  40.60
## 3   Sydney 151.21 -33.87

We add a "point geom" layer to our map:

ggplot(some_cities) +                                # pass the cities data to plot
    geom_sf(data = worldmap_data) +                  # layer 1: world countries
    geom_point(aes(x = lng, y = lat), color = 'red') # layer 2: points at city coord.

You may also add text labels for the cities. ggplot2 provides geom_text and geom_label (draws box around text):

ggplot(some_cities) +                                  # pass the cities data to plot
    geom_sf(data = worldmap_data) +                    # layer 1: world countries
    geom_point(aes(x = lng, y = lat), color = 'red') + # layer 2: points
    geom_label(aes(x = lng, y = lat, label = name),    # layer 3: labels
               hjust = 0, vjust = 1, nudge_x = 3) +    # labels appear next to point
    theme(axis.title = element_blank(),
          axis.text = element_blank())                 # disable the axis labels

Suppose you have observations that refer to geographic locations, e.g. countries, cities, neighborhoods. You may want to:

• link that data with existing information on that geographic entity e.g. country GDP, city population, poverty rate in the neighboorhood
• visualize that data geographically

For both things you need to:

1. agree on geographic (administrative) level for matching (e.g. country, state, municipality)
2. find identifiers to match the observations w/ respective geographic entities

TODO: schematische darstellung

• different names may refer to the same country (e.g. Germany / Federal Republic of Germany, BRD, etc.) → often not a good identifier
• ISO 3166-1 designates every country a two- or three-letter code (e.g. DE / DEU)
• often used in datasets (e.g. from the UN)

The Nomenclature of Territorial Units for Statistics (NUTS) divides the EU territory into regions at 3 different levels for socio-economic analyses of the regions.

Candidate countries and countries inside EFTA also have NUTS regions (Norway, Switzerland, Albania, Turkey, etc.).

In Germany, regions are hierarchically structured and identified by "Amtlicher Gemeindeschlüssel" (AGS ~ "municipality identificator"). The Gemeindeverzeichnis contains all regions with their name, AGS, area, population and other features.

Since 2006, Berlin uses "Lebensweltlich orientierte Räume" (LOR) to structure the city area into sub-regions at different levels, each of them identified by "LOR" code.

1. Prognoseräume (60 regions)
2. Bezirksregionen (138 regions)
3. Planungsräume (447 regions)

There are numerous file formats that store data for geographic information system (GIS) software. The most popular include:

• ESRI shapefile ("SHP"):
  • consists of at least three files in the same directory with filename extensions .shp, .shx and .dbf
• GeoJSON (.geojson, .json) and TopoJSON (.topojson, .json)
• KML (.kml, .kmz): Keyhole Markup Language, used by Google Earth
• GML (.gml): Geography Markup Language
• AutoCAD DXF (.dxf)

You may also encounter other formats such as vector data files (.svg, .poly) or databases (.sql, .sqlite).

Good news: sf can handle all popular file formats.

Two datasets:

• to combine these datasets, the "city" column might be used as identifier for matching
• not all cities in "pm" also appear in "city_coords" and vice versa

Datasets can also be joined with merge, but I find dplyr easier to use.

left_join(a, b, by = <criterion>): always retains rows on the "left side" and fills up non-matching rows with NAs.

library(dplyr)
left_join(pm, city_coords, by = 'city')

##        city pm_mg      lng      lat
## 1     Amman   999 31.91034 35.94846
## 2  Saltillo   869  1.00000  1.00000
## 3      Usak   814  1.00000  1.00000
## 4 The Bronx   284       NA       NA
## 5     Seoul   129  1.00000  1.00000

right_join(a, b, by = <criterion>): always retains rows on the "right side" and fills up non-matching rows with NAs. How many rows do you expect for a right join between pm and city_coords? How many of them will contain an NA value?

right_join(pm, city_coords, by = 'city')

##       city pm_mg      lng      lat
## 1    Amman   999 31.91034 35.94846
## 2 Saltillo   869  1.00000  1.00000
## 3   Berlin    NA  1.00000  1.00000
## 4     Usak   814  1.00000  1.00000
## 5 New York    NA  1.00000  1.00000
## 6    Seoul   129  1.00000  1.00000

inner_join(a, b, by = <criterion>): only retains rows that match on both sides.

How many rows do you expect for an inner join between pm and city_coords?

inner_join(pm, city_coords, by = 'city')

##       city pm_mg      lng      lat
## 1    Amman   999 31.91034 35.94846
## 2 Saltillo   869  1.00000  1.00000
## 3     Usak   814  1.00000  1.00000
## 4    Seoul   129  1.00000  1.00000

In a choropleth map, regions are shaded according to the measurement of some variable, e.g.:

• data downloaded from the Berlin Senate Dept. for Urban Dev. and Housing via FIS Broker
• includes indicators for monitoring of social development ("Monitoring Soziale Stadtentwicklung 2017")¹:

bln_sozind <- read_sf('excercises/data/bln_plr_sozind.geojson')
head(bln_sozind)

## Simple feature collection with 6 features and 10 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 386668.2 ymin: 5817761 xmax: 390764.3 ymax: 5820432
## epsg (SRID):    25833
## proj4string:    +proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## # A tibble: 6 x 11
##   PLANNAME EW2016 STATUS1 STATUS2 STATUS3 STATUS4 DYNAMO1 DYNAMO2 DYNAMO3
##   <chr>    <chr>    <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>
## 1 Stülers… 3114      4.98    1.52    8.16   20.5    -0.19    0.14   -0.42
## 2 Großer … 192      NA      NA      NA      NA      NA      NA      NA   
## 3 Lützows… 5518      6.83    2.08   17.7    34.3    -0.5    -0.03   -0.8 
## 4 Körners… 4704      6.27    1.49   20.7    45.9    -1.33   -0.76   -1.66
## 5 Nördlic… 1135      2.06    0       1.5     3.55   -0.77   -0.74   -0.42
## 6 Wilhelm… 2190      3.64    0.89    8.45   27.3    -1.23   -0.59    0.13
## # … with 2 more variables: DYNAMO4 <dbl>, geometry <MULTIPOLYGON [m]>

Directly plotting the map with fill color depending on STATUS4:

ggplot() + geom_sf(aes(fill = STATUS4), data = bln_sozind)

Problem: STATUS4 is continuous and so is the color scale. For choropleth maps, discrete ranges are better for discerning the colors.

Put the percentages into 5 discrete bins, change the color palette, adjust the appearance:

bln_sozind$child_pov_bins <- cut(bln_sozind$STATUS4, seq(0, 100, by = 20))
ggplot() + geom_sf(aes(fill = child_pov_bins), data = bln_sozind) +
    scale_fill_brewer(palette = 'OrRd', na.value = "grey90",
                      guide = guide_legend(title = 'Child pov.\nranges in pct.')) +
    coord_sf(datum = NA) +  # disable lines ("graticules") in the background
    theme_bw() + theme(axis.text = element_blank(), axis.title = element_blank(),
                       axis.ticks = element_blank())

Most of the time, you have at least two datasets: one containing the measurement(s) you want to show, the other containing the geo-data.

Example dataset: People at risk of poverty or social exclusion by NUTS 2 regions (tgs00107) from Eurostat.

pov_risk <- read.csv('excercises/data/tgs00107_pov_risk_nuts2.csv',
                     stringsAsFactors = FALSE)
pov_risk$risk_pct_bins <- cut(pov_risk$risk_pct, seq(0, 100, by = 10))

pov_risk_2016 <- filter(pov_risk, year == 2016)   # 2016 has fewest NAs
head(pov_risk_2016)

##   nuts year risk_pct risk_pct_bins
## 1   AT 2016     18.0       (10,20]
## 2 AT11 2016     15.1       (10,20]
## 3 AT12 2016     13.0       (10,20]
## 4 AT13 2016     26.0       (20,30]
## 5 AT21 2016     16.3       (10,20]
## 6 AT22 2016     18.4       (10,20]

We load the GeoJSON dataset for NUTS level-2 regions provided by Eurostat:

nutsrg2 <- read_sf('excercises/data/nutsrg_2.json')
st_crs(nutsrg2) <- 3857  # set the correct CRS
head(nutsrg2)

## Simple feature collection with 6 features and 2 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 1347007 ymin: 4105701 xmax: 3847132 ymax: 6631069
## epsg (SRID):    3857
## proj4string:    +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs
## # A tibble: 6 x 3
##   id    na                                                         geometry
##   <chr> <chr>                                            <MULTIPOLYGON [m]>
## 1 CY00  Κύπρος      (((3592720 4172836, 3594320 4178301, 3607520 4169324, …
## 2 CZ01  Praha       (((1588619 6463231, 1592220 6468696, 1617821 6476893, …
## 3 CZ02  Střední Če… (((1686224 6537392, 1685024 6526853, 1686624 6506167, …
## 4 CZ03  Jihozápad   (((1492615 6461670, 1536217 6432006, 1539017 6395706, …
## 5 CZ04  Severozápad (((1612621 6534270, 1600220 6515144, 1520216 6485089, …
## 6 CZ05  Severových… (((1729426 6582279, 1758628 6576034, 1786629 6556127, …

Both datasets contain a NUTS level-2 identifier so we can join them:

pov_risk_2016_geo <- left_join(nutsrg2, pov_risk_2016, by = c('id' = 'nuts'))
head(pov_risk_2016_geo)

## Simple feature collection with 6 features and 5 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 1347007 ymin: 4105701 xmax: 3847132 ymax: 6631069
## epsg (SRID):    3857
## proj4string:    +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs
## # A tibble: 6 x 6
##   id    na      year risk_pct risk_pct_bins                        geometry
##   <chr> <chr>  <int>    <dbl> <fct>                      <MULTIPOLYGON [m]>
## 1 CY00  Κύπρος    NA     NA   <NA>          (((3592720 4172836, 3594320 41…
## 2 CZ01  Praha   2016     10.1 (10,20]       (((1588619 6463231, 1592220 64…
## 3 CZ02  Střed…  2016     10.3 (10,20]       (((1686224 6537392, 1685024 65…
## 4 CZ03  Jihoz…  2016     10   (0,10]        (((1492615 6461670, 1536217 64…
## 5 CZ04  Sever…  2016     19.5 (10,20]       (((1612621 6534270, 1600220 65…
## 6 CZ05  Sever…  2016     11.7 (10,20]       (((1729426 6582279, 1758628 65…

(eu_pov_plt <- ggplot() + geom_sf(aes(fill = risk_pct_bins), data = pov_risk_2016_geo) +
    scale_fill_brewer(palette = 'OrRd', na.value = "grey90",
                      guide = guide_legend(title = 'Pct. of people at\nrisk of poverty')) +
    labs(title = 'Poverty in the EU 2016',
         subtitle = 'Percentage of people at risk of poverty or social exclusion',
         caption = 'source: Eurostat') +
    coord_sf(datum = NA) +  # disable lines ("graticules") in the background
    theme_bw() + theme(axis.text = element_blank(), axis.title = element_blank(),
                       axis.ticks = element_blank()))